arXiv:1505.07529v3 [math.NA] 6 Apr 2017 


Gaussian-Like Immersed-Boundary Kernels with Three Continuous Derivatives 

and Improved Translational Invariance^ 


Yuanxun Bao^’*, Alexander D. Kaiser^, Jason Kaye^, Charles S. Peskin^ 

“Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY, 10012, USA 


Abstract 

The immersed boundary (IB) method is a general mathematical framework for studying problems involving 
fluid-structure interactions in which an elastic structure is immersed in a viscous incompressible fluid. In the 
IB formulation, the fluid described by Eulerian variables is coupled with the immersed structure described 
by Lagrangian variables via the use of the Dirac delta function. From a numerical standpoint, the Lagrangian 
force spreading and the Eulerian velocity interpolation are carried out by a regularized, compactly supported 
discrete delta function, which is assumed to be a tensor product of a single-variable immersed-boundary 
kernel. IB kernels are derived from a set of postulates designed to achieve approximate grid translational 
invariance, interpolation accuracy and computational efficiency. In this note, we present new 5-point and 
6-point immersed-boundary kernels that are and yield a substantially improved translational invariance 
compared to other common IB kernels. 

Keywords: Immersed boundary method, fluid-structure interaction, discrete delta function, 
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1. Introduction 

The immersed boundary (IB) method was proposed to study flow patterns around heart valves O. In 
the IB formulation, a viscous incompressible fluid described by Eulerian variables is assumed to occupy 
the entire domain, which contains an immersed structure, described by Eagrangian variables, that moves 
with the fluid and exerts a force on the fluid. In the spatially discretized setting, the fluid domain is repre¬ 
sented by a uniform Eulerian grid and the immersed structure is configured as a collection of Eagrangian 
points or markers. The IB kernel plays a key role in communicating between the Eulerian and Eagrangian 
grids by spreading applied forces to the fluid and interpolating Eagrangian marker velocity. There are three 
main criteria for constructing an ideal IB kernel: grid translational invariance, interpolation accuracy and 
computational efficiency. It is a desirable property of an IB kernel to perform force spreading and velocity 
interpolation that are independent of the position of Eagrangian markers relative to the Eulerian computa¬ 
tional grid. In this case, if the IB method were applied to a translation-invariant linear system like the Stokes 
equations on a periodic domain, the results would remain the same despite shifts in position of Eagrangian 
markers relative to the Eulerian grid lJ3i|. There are functions that might serve as candidates for an IB kernel 


*This manuscript is an update to the the published article jT) with a new 5-point IB kernel. 

’Corresponding author 

Email addresses: billbaoOcims. nyu. edu (Yuanxun Bao), kaiserOcims. nyu. edu (Alexander D. Kaiser), 
jkayeOcims .nyu. edu (Jason Kaye), peskinScims .nyu. edu (Charles S. Peskin) 


1 





in terms of exact grid translational invariance; for example, the sine function sin(x)/x. However, the sine 
function is not computationally efficient because its support is unbounded. In fact, as we discuss later, exact 
grid translational invariance is inconsistent with the assumption of compact support [J]. The process of 
constructing a computationally efficient IB kernel that simultaneously has good interpolation accuracy and 
translational invariance is non-trivial. 

In the IB method, the 3D discrete delta function is assumed to be represented by a tensor product of a 
single-variable kernel (p{r). 



( 1 . 1 ) 


where x\,X 2 ,x^ are the Cartesian components of x and h is the meshwidth. This representation is not 
essential, but it significantly simplifies fhe discussion, since fhe single-variable kernel ^(r) is fhe objecf of 
inferesf. We firsf require fhaf 0(r) be confinuous for all real r and have compacf supporf, i.e., (p{r) - 0 for 
\r\ > Vs, where is fhe radius of supporf. Confinuify of (p is assumed in order fo avoid sudden jumps in fhe 
inferpolafed velocify or applied force as Lagrangian markers move fhrough fhe Eulerian grid. If furns ouf 
fhaf mosf IB kernels are even fhough fhe higher regularify is nof explicifly assumed. The reason for fhaf 
is sfill a mysfery, buf higher regularify of fhe IB kernel is a nice feafure fo have in certain applications, such 
as fhe inferpolafion of derivafives or fhe spreading of a force dipole. Compacf supporf of (p is required for 
compufafional efficiency. 

The function (p{r) is consfrucfed by requiring a subsef of fhe following momenf conditions: 


(i) 

Zerofh momenf: 

^4>{r- j) = 1, 

j 

(ii) 

Even-odd: 

2] (l>{r -j)= Yj ~ 

j even / odd 

(iii) 

Eirsf momenf: 

- j) (p{r - j) ^ 0, 

j 

(iv) 

Second momenf: 

- j)^(p{r - j) = K, for some consfanf K, 

j 

(V) 

Third momenf: 

- if(p{r - i) ^ 0. 


,/■ 


The mofivafion of imposing momenf conditions is well discussed in |l3l|4l. Briefly, fhe zerofh momenf con¬ 
dition implies fhaf fhe fofal force is fhe same in fhe Eulerian and Eagrangian grids when 6h is used for force 
spreading. The even-odd condifion implies (i), and was originally proposed fo avoid fhe “checkerboard” in- 
sfabilify fhaf may arise from using a collocafed-grid fluid solver. Eiu and Mori 14;] generalized fhis condifion 
fo fhe so called “smoofhing order” condifion and showed fhaf if has fhe effecf of suppressing high-frequency 
errors and preventing Gibbs-fype phenomena. Conservafion of fofal forque reties on fhe firsf momenf con¬ 
difion. Moreover, (i) and (iii) guaranfee fhaf a smoofh funcfion is inferpolafed wifh second-order accuracy 
when 6h is used for inferpolafion. The second momenf condifion wifh A' = 0 and fhe fhird momenf condifion 
are needed fo derive kernels wifh a higher order of inferpolafion accuracy. 

In addifion fo momenf conditions, ^{r) is required fo safisfy fhe sum-of-squares condifion, 

^ {(p{r - j))^ = C, for some consfanf C. (1.2) 
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The sum-of-squares condition Eq. (1.2)| is a weaker version of exact grid translational invariance, 

^(ri, r 2 ) ^ ^ 0(ri - j) (p{r 2 - j) = 0(ri - ^ 2 ), for all n, r 2 . (1.3) 


In other words, the coupling of 0(r) between any arbitrary two points r \, r 2 is a function of ri - r 2 only. 
However, it can be shown that |Eq. (1.3) is inconsistent with the assumption of 0 being compactly supported 
ll3l . The sum-of-squares condition does give some information about the coupling of cf), since it can be 
deduced from the Cauchy-Schwarz inequality that 


\Kn,r2)\ 


“^(^2 - i) 

j 


< C, for all ri,^ 2 . 


(1.4) 


Eq. (1.4) guarantees that the coupling between two Eagrangian markers is strongest when the markers 


coincide, and furthermore Eq. (1.2)|implies that the self-coupling is independent of the marker position. 


IB Kernel 

Support 

rs 

Even-Odd 

Zeroth 

Moment 

Eirst 

Moment 

Second 

Moment 

Third 

Moment 

Sum of 
Squares 

Regularity 

Standard 

3-point 

1.5 

X 

/ 

/ 

X 

X 

1 

2 


Smoothed 

3-point 

2 

X 

/ 

/ 

X 

X 

X 


Standard 

4-point 

2 

/ 

/ 

/ 

X 

X 

3 

8 


Smoothed 

4-point 

2.5 

X 

/ 

/ 

X 

X 

X 


Standard 

6-point 

3 

/ 

/ 

/ 

0 

/ 

67 

128 


New 

5-point 

2.5 

X 

/ 

/ 

38 

60 60 

/ 

« .393 


New 

6-point 

3 

/ 

/ 

/ 

59 V 29 

60 20 

/ 

» .326 

^3 


Table 1: Common immersed-boundary kernels with their properties and moment conditions they satisfy. /: the kernel satisfies 
the moment condition; X: the kernel does not satisfy the moment condition. In the second moment column, the value of the 
second moment constant K is given when the second moment condition is satisfied. The regularity column shows the number of 
continuous derivatives each IB kernel has. 


In Table 1 we list some common IB kernels and the conditions they satisfy. The most widely used IB 
kernel is the standard 4-point kernel f3i|. The standard 3-point kernel satisfies the zeroth moment condition 
but not the even-odd condition. It was first introduced in an adaptive IB method using the staggered-grid 
discretization iSl. The standard 6-point kernel (with K = 0) satisfies all the moment conditions listed above 
Il6l . It can be shown that the standard 6-point kernel interpolates cubic functions exactly and smooth func¬ 
tions with fourth-order accuracy. However, as shown in Eig. 4[ it has a larger deviation from translational 
invariance for a pair of markers with distance d w 2.5 even compared to the standard 4-point kernel. In 
terms of its defining postulates, our new 6-point kernel differs from the standard 6-point kernel only in the 
nonzero second-moment constant K (the sum-of-squares constant C is determined once K is fixed). The new 
5-point kernel assumes the same postulates as the new 6-point kernel except for the “even-odd” condition. 
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The special choice of K given in Eq. (2.12) and Eq. (2.15)] lead to new 5-point and 6 -point kernels that are 
and significantly improve translational invariance compared with other IB kernels. The construction of 
an IB kernel with a positive and constant second moment K was originally motivated by the important phys¬ 
ical implication of the second moment in particle suspensions, namely it is associated with the quadrupole 
correction in the Eaxen relation for a rigid sphere in an arbitrary Stokes flow fTl- The result that the new 
kernels have three continuous derivatives is unexpected, however, this makes the new kernel more generally 
useful. By applying a smoothing technique to the standard IB kernels, Yang, et al. [8 ] developed a family 
of IB kernels whose first derivative satisfies up to the second moment condition for the derivative. They 
showed that these derivative moment conditions are intrinsically linked to the error of force spreading in the 
IB scheme, and IB kernels that satisfy these conditions can significantly reduce non-physical spurious oscil¬ 
lations of force spreading in moving-boundary problems. By differentiating the moment conditions satisfied 
by ^{r), we can verify that the derivative of our new kernels satisfy up to third moment conditions as 
advocated by Yang et al. [S ]. We will also include the smoothed 3-point and 4-point kernels in the compar¬ 
ison of translational invariance in section 3 Liu and Mori developed a MATLAB routine that automatically 
generates all the standard IB kernels as well as many others |5l- We have also made our MATLAB codes 
for generating the new kernels available at https: //github. com/stochasticHydroTools/IBMethod 


2. Two new kernels 

2.1. A new 5-point kernel 


Our new 5-point kernel satisfies the sum-of-squares condition |Eq. (1.2) and the moment conditions (i), 
(iii)-(v), but it does not satisfy the even-odd condition (ii). The support is defined to be five grid points, i.e.. 


rs = ^- We follow a similar derivation of the standard 4-point kernel [T]. By first restricting re [- 2 ^ 
have 5 unknowns: 


we 


{(pir - 2 ), (f>{r - 1 ), 0 (r), (f>{r ■+■ 1 ), 0 (r -F 2 )). 


(2.1) 


Note that the moment conditions (i), (iii)-(v) are four linear equations in these five unknowns, and we can 
express all the other four unknowns in terms of (p{r), 


(p{r “ “ Y 2 + r^ + 2r^ - r - I'j, 

(f){r - 1) = 7 (-4(l>{r) - 3Kr - K - r^ - r^ 4r + , 

6 ' ' 

(f){r + 1) = 7 {-4(p{r) + 3Kr - K r^ - r^ - 4r 4), 
6 ^ ' 

(P{r + 2) = (20(r) - 3Kr + 2K - ? + 2r^ + r-l). 


( 2 . 2 ) 

(2.3) 

(2.4) 

(2.5) 


Eor the special value r = 1/2 so that <;& = 0, we get 


1 


1 


1 


1 


1 


n-7;\, n7:\, nTtW = {t7^4K - - 4K),—{9 - 4K),—{4K -\) \. 


16 


16 


16 


16 


( 2 . 6 ) 


Substituting these values into the sum of squares condition Eq. (1.2)[ we obtain an expression for C{K), 

(2.7) 
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where the value of K remains to be determined. Next we solve the quadratic equation of ^(r) from the 
sum-of-squares condition Eq. (1.2) by using Eqs. (2.2)1to|(2.5)|and|Eq. (2.7)[ 


m = ^ (-40/: - 40r^ + 136 + ^j2/3{r) + 2y{r )) , 


( 2 . 8 ) 


where 




m = 

y(r) = -40r2(35/-202r2 + 31l). 


S400K/ + 256S0Kr'^ 


6840/:+ 3123, 


(2.9) 

( 2 . 10 ) 


We note that the positive square root is chosen in|Eq. (2.^ because of th e continui ty assumption of 0, i.e., 


by setting r = we select the branch that gives (p = ^{9 - 4K) as in Eq. (2.6) 

We have the freedom to choose K to construct 0 with higher regularity. A symbolic calculation in 
Mathematica matching derivatives of (f>{r) at r - ^ reveals that for any K e |^0, 0 e Eor the second 

derivative to be continuous at r = I, K must satisfy 


720/:^ -912/:+ 275 - 0. 


( 2 . 11 ) 


We can verify (by plotting) that exactly one of the two roots of Eq. (2.11) guarantees that f3{r) + 7 (r) > 0, 
so that (p{r) is real-valued for r e [ 5 , jj, and this special value of K is 


^=^(38-V69). 


( 2 . 12 ) 


If we proceed further with matching the third derivative to be continuous at r = i, then K must satisfy 


60480/:^ - 116208/:2 + 73260/: - 15125 = 0, 


whose roots are 


K = 


M' 60 


1(38+V«), 1(38- 


(2.13) 


(2.14) 


We observe that our choice of K in Eq. (2.12) is among the roots of the third derivative matching condition, 


and therefore, it also makes the third derivative of cp continuous at r = 5 . 

Remarkably, the derivative matching conditions at r = ^ are sufficient to ensure that cp everywhere. 
Existence and continuity of derivatives was never assumed a priori, but was only a consequence of an 
appropriate choice of K. Moreover, note that for /: e |^0, ^ ^38 - we have (f) < 0. Since 

^(i) “ ( 2 ) “ this implies that (p{f) takes negative values in a neighborhood of r = 5 . We emphasize 

that the special choice of K given by Eq. (2.12) is the smallest positive K for which (p is non-negative, and 
it is also the only value of K that gives three continuous derivatives of 0. 

2.2. A new 6-point kernel 


Our new 6 -point kernel satisfies the sum-of-squares condition Eq. (1.2) and the moment conditions 
(ii)-(v) (and therefore (i)) with the second-moment constant 


K = 


59 


(2.15) 




















The derivation of this kernel follows a nearly identical procedure to that of the new 5-point kernel. First, 
the sum-of-squares constant C can be expressed in terms of K by considering the special case r = 0. Next, 
by restricting r to the interval [0,1], we have 6 unknowns: - 3), 0(r - 2), (t>{r - 1), (p{r), (p{r + 1), (j){r + 2) 

and 6 equations (the even-odd condition accounts for two equations). By expressing all the other unknowns 
in terms of ^{r - 3) using (ii)-(v), we can solve for - 3) from the quadratic equation determined by 
Eq. (1.2)1 The continuity assumption of (p is now used to select the appropriate root to piece together a 
continuous function, i.e., by setting r = 0, we select the root that gives (/'(-3) = 0. As mentioned earlier, (p 
being follows implicitly from our defining postulates, i.e., (p'{-3) = 0. We have the freedom to choose 
K so that (p"{-3) = 0, which uniquely determines the special value Eq. (2.15) The formula for the new 
6 -point kernel is given by 


Q 99 

/3(r) =---{K^r^)r + {- 


7 T 

■lK)r--r^ 


1 

T2 


{i3K - l)r + Pf + ((4 - 3K)r - r^) 


(2.16) 


(2.17) 


(p{r-3) ^ 


-P{r) -I- sgn (I - is:) 


56 


1 K + {3K- 

4>{r-2) ^ + + ^ 


l)r 


1 (4 

<p{r - 1) ^ l(p{r - 3) -H - — 


3K)r 


<p(r) = 2<p(r 
(p{r -I- 1) = - 3(p{r 


3) + ? 


6 

K + r'^ 


r 

~6 




4 

(4 


(p{r + 2) - (p{r ■ 


1 K + r^ 
”■16 ■"“ 8 “ 


- 3K)r 
6 6 
(3K - l)r 


12 


r-^ 

12 


(2.18) 

(2.19) 

( 2 . 20 ) 
( 2 . 21 ) 
( 2 . 22 ) 
(2.23) 


Note that, in the formula presented above, r € [0,1]. The new 5-point and 6-point kernels are Gaussian-like 


function, as shown in Eig. 1 and Eig. 2 and they both have three continuous derivatives. As a comparison. 


the standard 3-point, 4-point, 6-point kernels and their continuous first derivative are plotted in Eig. 3 We 
notice that the new 5-point and 6-point kernels are non-negative for all r, whereas the standard 6-point 
kernel has negative tails. 


3. Numerical Tests 


In this section, we demonstrate a significant improvement in the translational invarianc^ of our new 
6 -point kernel. We randomly select 10^ pairs of Eagrangian markers Xi, X 2 in a periodic box [0,32]^ with 


meshwidth h = \ and compute the 3D generalization of Eq. (1.3) 


d(Xi,X 2 ) = ^ 6h{x - Xi) Shi'S. - X 2 ), 


(3.1) 


xegii 


*The test that we use actually checks for rotational invariance at the same time, since it involves the Euclidean distance between 
a pair of markers, and not merely the vector from one marker to the other. 


6 





















-3-2-10 1 2 3 


r 

(a) 



r 

(b) 


Fig. 1: (a) The new 5-point kernel compared with the Gaussian 
derivatives of the new 5-point kernel. 


with the second moment given by Eq. (2.12) 


(b) The first three 
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Fig. 2: (a) The new 6-point kernel compared with the Gaussian 
derivatives of the new 6-point kernel. 


with the second moment given by Eq. (2.15) 


(b) The first three 
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Fig. 3: (a) The standard 3-point, 4-point, and 6-point kernels, (b) The first derivatives of the standard 3-point, 4-point, and 6-point 
kernels. 


where x denotes a grid point on the Eulerian grid gh. In Fig. 4 we plot d(Xi, X 2 ) normalized by the constant 
from Eq. (1.2) versus the distance d = |Xi - X 2 I. The data are binned according to d - |Xi - X 2 I, 
and error-bars showing the maximum, mean and minimum of each bin are overlaid with the data. If an IB 
kernel were exactly translation-invariant, the plot of d(Xi, X 2 ) would be a curve. The spreading pattern in 
the data around this curve clearly indicates that none of the IB kernels compared here are exactly translation- 
invariant, but gives a qualitative picture of how close to translational invariance each kernel is. The data 
of the new IB kernels and the smoothed 4-point kernel almost form a curve, while the data of the other 
kernels have larger deviations from the mean. Moreover, the deviation in the data of the new IB kernels is 
uniform for all distances, but the standard 6-point kernel has a much larger deviation for <i w 2.5. For a more 
quantitative comparison, we summarize the maximum standard deviation of all bins for each IB kernel in 


Table 2 The maximum standard deviation of the new IB kernels is an order of magnitude smaller than that 


of the standard IB kernels, and is about half of the deviation of the smoothed 4-point kernel. 

In terms of computational cost, we summarize the computation time of §(Xi,X 2 ) for 10^ pairs of Ea- 
grangian markers using the kernels we have compared. The timings are based on simulations performed on 
a desktop with Intel Core 17-4770 CPU 3.40GHz under the MATEAB environment. The main cost of using 
an IB kernel in spreading/interpolation depends on its support size. In Table 2 the new 5-point kernel is 
about two times more expensive, and the new 6-point kernel is about three-to-four times more expensive 
than a 4-point kernel in our comparison, because the new 5-point and 6-point kernels communicate with 
125 and 216 nearby grid points respectively in spreading/interpolation in 3D, while a 4-point kernel only 
communicates with 64 nearby grid points. The smoothed 3-point and 4-point kernels are more expensive 
than their standard counterparts in that they have wider supports as shown in Table 1 In all respects, the 


new IB kernels achieve significant improvement in grid translational invariance with a modest increase in 
computational cost. 


4. Conclusion 

In this note, we have presented new immersed-boundary kernels that are used for force spreading and 
velocity interpolation in the immersed boundary method. The new IB kernels are distinguished from other 
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d= |Xi -Xal 

(c) 


Fig. 4: Normalized 5(Xi,X2) is plotted versus d = |Xi - X 2 I for 10^ pairs of randomly selected Lagrangian markers. The data are 
binned according tod = |Xi -X 2 I and error-bars showing the maximum, mean and minimum of each bin are overlaid with the data. 
The deviation in the data gives a quantitative measure of translational invariance of an IB kernel. The standard deviation in the 
data for the new 6-point kernel (blue) is an order of magnitude smaller than for the standard IB kernels, (a) The standard 6-point 
kernel vs. the new 5-point kernel vs. the new 6-point kernel, (b) The standard 3-point kernel vs. the standard 4-point kernel, (c) 
The smoothed 3-point kernel vs. the smoothed 4-point kernel. 
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Standard 

3-point 

Smoothed 

3-point 

Standard 

4-point 

Smoothed 

4-point 

Standard 

6 -point 

New 

5-point 

New 

6 -point 

maximum 

std. dev. 

0.0428 

0.0212 

0.0168 

0.0083 

0.0296 

0.0051 

0.0042 

computation 

6.34s 

9.52s 

9.01s 

12.06s 

31.19s 

17.54s 

30.86s 


time 


Table 2: Maximum standard deviation of d(X|,X 2 ) over all bins for various IB kernels, and the computation time for computing 
d(Xi,X 2 ) for 10^ markers. 


existing IB kernels by their nonzero second-moment constant K. A special choice of K leads to a 5-point 
or 6-point IB kernel that is and features substantially improved translational invariance compared with 
the existing standard IB kernels. Recently, we have successfully applied the new IB kernels to a new IB 
method with an exactly divergence-free interpolated velocity field [9|, in which derivatives of the discrete 
delta function are involved, and have achieved 10^ to 10^ times improvements in volume conservation of 
the IB method. We believe that, in many other applications in which derivatives of the IB kernel are needed, 
the improved grid invariance and regularity of the new IB kernels will be worth its extra computational cost. 
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